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1.  Introduction 

Qubit  unitary  algorithms  are  a  novel  computational  scheme  for  the  solution  of  nonlinear 
problems  of  physics.  Not  only  are  they  idealized  parallelized  on  classical  supercomputers  but  will 
be  immediately  applicable  to  quantum  computers  when  they  come  on  line. 

While  a  binary  bit  can  take  on  only  the  values  “0”  or  “1”,  the  qubit  resides  on  a  surface  of 
the  unit  sphere  with  state  |<y)=a0|0)  +  1 1) ,  a  superposition  of  states  |0)  and  |l) .  The 

(complex_  probability  amplitudes  |a0|2  +|aj|2  =  1 . 


Bit 

0 


1 


Pbit 

0 


Qubit 

0 
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Fig.  1  A  schematic  of  a  bit,  a  probability-bit  ‘Pbit’,  and  a  qubit. 


Things  get  more  interesting  when  one  considers  the  interaction  of  two  qubits : 


kl^2)  =  «Ool00>  +  «Oll01>  +  alol10)  +  «lll11>  -  With  S 


aU 


=  1 


Of  much  importance  in  quantum  information  is  qubit  correlations  -  i.e.,  quantum  entanglement.  It 
is  this  feature  that  permits  quantum  information  to  solve  problems  like  n-digit  factorization  into  two 
prime  number  that  scale  polynomial^  in  n  (Shor  algorithm)  while  on  classical  supercomputers  the 
scaling  is  exponential  in  n. 


The  typical  quantum  lattice  algorithm  for  the  evolution  of  a  scalar  mesoscopic  field  is  the 
introduction  of  2  qubits  per  lattice  node.  One  then  introduces  an  appropriate  sequence  of  unitary 
collide-stream  operations  that  entangle  the  on-site  qubits  and  then  streams  this  entangled 
information  to  neighboring  lattice  sites.  This  sequence  of  collide-stream  unitary  operations  has 
been  mimicked  in  some  successful  NMR  experiments  by  the  Cory  group  [xxx].  In  particular,  since 
a  unitary  operator  U  is  related  to  a  Hamiltonian  H  by  U  =  exp (iH) ,  the  NMR  experimentalists  and 

readily  determine  the  NMR  pulses  and  their  durations  in  order  to  emulate  the  desired  unitary 
operations.  Thus,  from  the  qubit  unitary  algorithm  for  the  one-dimensional  (ID)  Burgers  equation 

du  du  d2u 

— - 1"  U  ~ —  —  V - FT 

dt  dx  dx 

one  can  immediate  determine  the  required  electromagnetic  pulses  and  polarizations.  In  2004, 
using  16  lattice  sites,  with  2  qubits/site,  the  Cory  group  could  reproduce  the  early  stage  of  the 
steepening  of  a  sinusoidal  wave  front  towards  a  shock.  The  dashed  points  are  experimental  data 
while  the  continuous  curve  is  the  analytic  solution  to  the  Burgers  equation 
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Fig.  2  A  comparison  between  theoretical  (lines)  and  NMR  data  points  in  the  early  stages  of  shock 
development  in  the  solution  of  the  ID  Burgers  equation.  The  NMR  pulse  sequence  is  determined 
from  the  Flamiltonian  representation  of  the  unitary  collide-stream  operators  of  the  qubit  algorithm 
that  solves  the  Burgers  equation. 


It  is  also  of  interest  to  consider  some  of  the  ideas  being  put  forward  for  a  quantum  computer.  In 
particular  the  Zoller  group  (PRL  -  2008,  2011)  has  proposed  using  the  alkaline  earth  atom  87  Sr  in 
their  Bose-Einstein  condensate  optical  lattices. 


storage  lattice 


10)  |i) 


»  7 

mi  =  ~2  "»/  =  ”  j 


Fig.  3  Quantum  computing  with  independent  lattices:  (a)  qubits  in  long-lived  states  in  a  storage 
lattice  are  transferred  to  a  completely  independent  transport  lattice  for  operations  between  distant 
qubits,  or  addressed  individually  by  coupling  to  a  level  that  is  shifted  by  a  gradient  magnetic  field. 
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(b)  This  is  accomplished  by  encoding  qubits  in  nuclear  spin  states,  producing  independent  lattice 
for  the  lS0  and  3P0  levels,  and  using  3  P2  for  individual  addressing. 

One  can  readily  encode  qubits  in  the  nuclear  spin  states  with  different  magnetic  quantum  numbers 

ntj-  in  particular  one  can  define  |0}=  1S0,mI  =  -9^,  and  the  state  |l)=  1S0,mI  =-^, 

Fig.  3b.  These  are  quite  stable  states  since  they  are  very  insensitive  to  magnetic  field  fluctuations. 
Moreover  quantum  tunneling  is  negligible  resulting  in  deep  optical  lattices  for  storage  and  transport 
lattices.  We  thus  see  some  similarity  with  these  ideas  of  developing  a  quantum  computer  and  the 
use  of  unitary  operators  for  shifting  qubits  and  collisional  coupling  of  qubits. 


2.  ID  Quantum  Lattice  Algorithms  for  Classical  Solitons 

We  have  investigated  quantum  unitary  algorithms  for  both  the  ID  Korteweg-de-Vries  and 
the  Nonlinear  Schrodinger  equations  [G.  Vahala,  J.  Yepez  and  L.  Vahala,  Phys.  Lett.  A310,  187- 
196  (2003)].  In  particular  to  recover  the  ID  nonlinear  Schrodinger  equation  for  bright  solitons 

,d\i f  d2y/  |  ,2  ,,, 

=  v  (D 

we  introduced  2  qubits/lattice  site  and  employed  the  following  sequence  of  unitary  collide-stream 
operators  on  the  spinor 

M0}=f“)  (2) 

\PJ 

to  obtain  the  spinor  evolution  equation 

NLS:  \y/(t  +  At))=S2CS2C.S2CS2C .SfCS^  .SfCS^  \y/(t))  (3) 

where  the  unitary  5,  -operator  streams  the  probability  amplitude  a ,  while  S2  streams  /5 .  The 
unitary  collision  operator  C  is  known  as  the  square-root-of-swap: 


C  =  - 
2 


1 —i  1 +i 
1 +i  1 - i 


with  C=I 


(Here  /  is  the  identity  matrix). 


On  the  otherhand,  to  recover  the  KdV  equation 


dy/  dy/  d3y/ 

dt  dx  dx3 


=  0 


(6) 


with  the  nonlinear  convective  term  y/  dy/  /  dx ,  the  omission  of  the  standard  2nd  derivative  but  the 


appearance  now  of  the  3rd  derivative,  the  unitary  collide-stream  algorithm  interwines  the  streaming 
on  the  two  probability  amplitudes 


KdV:  \yf(t  +  &))=S]CT  ST2C  .STxCTS2C.SlCS2CT .SXCST2CT  \y/(t))  .  (7) 


The  unitary  collision  operator  is  now  the  square-root-of-NOT 
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KdV:  C  =  -j= 

V2 


1  1 

v  -i  1  y 


with  C=-I. 


(8) 


The  beauty  of  these  ID  solition  solutions  is  that  their  solutions,  though  nonlinear,  have  an 
arbitrary  scaling  parameter: 


NLS:  y/ NLS  (x ,t) — y[2  a  .Exp 


bx 

VT 


(  u  2 


-a 


.  Seek  a(y[2x-bt j 


(9) 


where  a  and  b  are  free  parameters,  controlling  the  soliton  amplitude  and  velocity.  For  the  KdV- 
soliton,  there  is  still  one  free  parameter  which  couples  the  soliton  amplitude  and  velocity. 


KdV:  y/Kdv(x,t)  =  ^a  Seek7 


\[a 


(. x-at ) 


(10) 


3.  Qubit  Unitary  Algorithms  for  3D  and  2D  Quantum  Turbulence 

Bose-Einstein  Condensates  (BECs)  of  weakly  interacting  dilute  gases  are  playing  an 
extremely  important  role  in  both  atomic  and  condensed  matter  research.  If  the  BEC  is  in  a 
magnetic  trap,  then  the  spin  is  aligned  to  the  magnetic  field  and  the  mean  field  representation  of 
the  evolution  of  the  ground  state  is  given  by  the  Gross-Pitaevskii  (GP)  equation  for  the  single-atom 
wavefunction 


.jfy '(*,*) 

dt 


^-V2-p+ g\y/(x,tf  y '(x,f) 


(11) 


p  is  the  chemical  potential  and  the  nonlinear  term  g\y/\2y/  arises  from  the  weak  s-wave  Boson 
interactions.  Under  the  Madelung  transformation,  the  GP  Eq.  (11)  can  be  written  in  fluid-like 
variables  of  mean  density  p(x,t)  and  mean  superfluid  velocity  us(x,r): 


y/(x,t)  =<Jp(xj)  exp[f<jp(x,r)] 
u  s{x,t)=V(p 


(12) 


with 


^  +  V-(pus)=0  ,  and  p 

at 


( 


\ 


^  +us-V 

at 


us  -  Vp  +  V  -Z 


(13) 


and  barotropic  pressure  :  p~gp2 .  Thus  the  GP-BEC  ground  state  corresponds  to  the 
dynamics  of  a  compressible ,  inviscid  Euler-like  superfluid,  but  with  a  novel  quantum  stress  tensor 


=  P 


d2\np 

dxj  dxj 


(14) 


It  was  finally  determined  that  the  Bose-Einstein  ground  state  wave  function  is,  in  fact,  an  order 
parameter,  with  a  topological  singularity  at  y/  =  0 .  This  topological  singularity  is  nothing  but  the 
core  of  a  quantum  vortex  with  mean  density  p  =  0 . 

Experimentally,  the  trapped  BEC  was  subject  to  rotation.  When  the  angular  momentum 
had  reached  a  threshold  value,  line  vortices  appeared  within  the  BEC  -  the  number  of  such 
vortices  being  dependent  on  the  value  of  the  angular  momentum  since  all  quantum  vortices  have 
the  same  strength 
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<£us  •  cU  =  me  (15) 

Indeed,  from  the  BEC  density  data,  one  finds  the  formation  of  a  vortex  lattice  as  the  BEC  is  rotated 
with  higher  angular  momentum: 


Fig.  4  The  formation  of  a  BEC  vortex  lattice  as  the  BEC  is  rotated  more  rapidly.  The  quantum 
vortex  can  be  readily  identified  by  the  zero  density  at  the  quantum  vortex  core,  which  in  these  plots 
of  the  BEC  density  show  up  as  cavities:  at  the  top,  the  BEC  rotation  is  too  low  for  the  existence  of  a 
quantum  vortex.  Moving  clockwise,  one  finds  16  quantum  cores,  then  70  quantum  cores  and 
finally  130  quantum  cores  [from  Ketterle  et.  at] 

There  are  considerable  differences  between  a  quantum  vortex  and  a  classical  vortex: 

(a)  a  classical  vortex  can  exist  in  an  incompressible  fluid,  but  the  quantum  vortex  is  essentially 
compressible  (since  p  =  0  at  the  vortex  core), 

(b)  the  presence  of  viscosity  is  essential  for  classical  vortex  reconnection,  while  quantum  vortex 
reconnection  occurs  without  viscosity  due  to  the  topological  singularity  at  the  vortex  core  and  the 
novel  quantum  stress  tensor  E. 

(c)  classical  vortices  (or  eddies)  can  have  a  continuous  range  of  values  for  circulation,  and  this 
play  a  significant  role  in  kinetic  energy  cascades:  in  3D,  energy  cascades  from  large  sale  eddies  to 
small  scale  eddies  till  viscosity  eventually  destroys  the  smallest  scale  eddies  and  dissipates  them 
into  heat.  On  the  otherhand,  quantum  vortices  are  identical,  with  the  same  quantized  circulation 
strength,  Eq.  (15).  Moreover,  it  is  energetically  unfavorable  for  the  persistence  of  high  winding 

number  vortices:  y/  =  y[p  exp(m<p)  ,n  >  2 .  These  degernate  winding-number  n  vortices  rapidly 

split  into  n  winding  number-1  vortices.  This  is  readily  seen  in  our  simulations,  Fig.  5,  which  show 
the  evolution  of  quantum  line  vortex  core  isosurfaces  for  a  degenerate  initial  winding  number 
n  =  5 .  Initially,  there  are  4  such  line  vortices  in  each  2D-plane  for  a  total  of  12  vortices,  Fig.  5a. 
The  bounds  on  the  color  scheme  on  the  vortex  core  isosurface  are:  blue  for  the  isosurface 
rp  =  0(mod  2k),  while  <p->2;r(mod  2k)  isosurface  is  in  red.  The  5-fold  periodicity  is 
evident.  On  the  walls,  we  plot  the  phase  information  of  the  wavefunction  as  it  intercepts  the 
boundaries.  From  this,  the  topological  singularity  (p  =  0)  acts  like  a  branch-point  singularity,  Eq. 

(12),  with  5  branch  cuts  connecting  each  quantum  vortex  core  intercept  with  the  boundaries.  Very 
rapidly,  Fig.  5b,  these  degenerate  isosurfaces  split  into  winding  number-1  vortices  (as  can  be 
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readily  seen  from  the  phase  information  on  the  walls  and  the  single  branch  cut  joining  the  singular 
vortex  cores.  The  oscillations  on  the  vortex  cores  are  due  to  the  excitation  of  Kelvin  waves  -  and 
this  we  will  discuss  later. 


(a)  t=  0  (b)  t  =  3000 


Fig.  5  The  splitting  of  the  5-fold  degenerate  quantum  vortices,  (a),  into  5  non-degenerate  quantum 
vortex  cores,  (b).  Lattice  grid  2048 3 .  Phase  information  on  both  the  topological  quantum  vortex 
cores  and  the  walls :  rp  =  0 (blue)  ,cp  =  In  [red).  The  oscillations  on  the  vortex  cores  as  seen  in 
(b)  are  due  to  Kelvin  waves. 

(d)  in  the  case  of  classical  vortices  the  evolution  equations  are  dissipative,  while  the  BEC 
quantum  vortex  evolution  is  determined  by  a  conservative  (non-dissipative)  Hamiltonian  system. 


Feynmann  put  forth  the  definition  that  quantum  turbulence  arises  with  a  3D  tangle  of 
quantized  vortices,  Fig.  6 


1 

plgt 

l&i. V-v  '*'>■ 

ml 

1  .  ^  - 
JiV  K .  -\x s  *  ^ V  :ttr  * 

I  ^  ^  V  \  fflh  '  f  v  .'  •  > v* :  ' 

1  (  *»  «<-  ’*  ^  .  V4 .  *7?  4. 

(a)  t=0  quantum  core  isosurfaces 


(b) 


later  turbulent  state  of  tangled  quantum  vortices 


Fig.  6  Evolution  of  quantum  turbulence,  defined  as  a  tangle  of  quantum  vortices. 
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4.  Qubit  Unitary  Lattice  Algorithm  for  Quantum  Turbulence 

To  develop  an  algorithm  for  the  3D  scalar  GP  Eq.  (11),  we  need  introduce  just  the 

minimal  2  qubits/spatial  node  with  its  Hilbert  space  of  dimension  22  -  i.e.,  the  states  |00),  |0l)  , 

1 10),  1 1 1)  -  and  focus  on  the  probability  amplitudes  for  the  two  states  1 10)  and  |01).  Let  us 
denote  these  probabilities  at  spatial  grid  point  x  at  time  t  as  a(x,t)  and  fi(x,t)  and  consider  the 
spinor  field 


MV) 


( a(x,r)  ' 
p(x,t)  , 


(16) 


We  now  consider  the  interleaved  unitary  collide-stream  operator  sequence  for  spatial  direction  x 


Ix,a  —  SfrxfjC  ,  (7—  0,1  ,  (17) 

where  the  unitary  collision  operators  that  entangles  locally  the  amplitudes  is  our  yfsWAP  -gate, 
Eq.  (4),  and  the  unitary  streaming  operator  that  spreads  this  entanglement  to  the  nearby  lattice 
sites  is  given  by 


Vo  Mx))= 


a(x  +  Axx) 

m  . 


ViMx))= 


a(x) 

P(x  +  Axx) 


For  3D,  the  full  unitary  evolution  operator  is 


c/[n]=4/,2,4exP 


—  l 


.  e2Q(x) 


j  2  r  2  r  2 

o4lexP 


—  l 


,e2Q(x) 


(18) 

(19) 


where  we  have  introduced  a  phase  rotation  factor  e2Q(x) ,  with  e  « 1 . 

The  mesoscopic  qubit  lattice  time  evolution  algorithm  is  thus 
\xi/(x,t  +  At))  =  t/[x,Q]  \y/{x,t)) 

In  the  spatial  continuum  limit,  Ax  =  0(e) ,  this  reduces  to  the  spinor  equation 

\x//(x,t  +  At))  =  |t^(x,r))  -z'e2[-t7xV2+Q(x)]  |yA(x,r))  +  o(e4) 
where  ax  is  the  Pauli  spin  matrix 


Thus,  under  diffusion  ordering, 

At-  Ax2  ~o(e2) 


(20) 

(21) 


(22) 


and  taking  moment  contraction  of  the  spinor  field  in  order  to  recover  the  scalar  order  parameter 
wavefunction  \//(x,t ): 

e->  1  (23) 

we  obtain  the  desired  GP  equation  for  y/(x,t) 


-  a  l-g|y/-(x,z)|2J  y/(x,t)  +  o(e2) 


dt 

on  identifying  Q  ^  g\y/\2 . 

•  Exact  Conservation  of  the  Mean  Density 


(24) 
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Since  the  qubit  lattice  algorithm  is  unitary,  the  norm  of  the  spinor  |i//)is  automatically 

conserved.  However,  there  can  be  fluctuations  in  the  mean  density  (p)  =  J  d3x\y/ (x,t)\2  if  there 
is  overlap  in  the  components  of  the  2-spinor 


(8p)  =  jd3x^\y/(x,tf  -(v'lv'))  =  jd3x(a*p  +aj8t) 


Thus,  if  a  is  purely  imaginary  and  /5  is  purely  real  (or  vice  versa),  then  the  overlap  between  the 
two  components  vanish  and  (Sp)  =  0 . 

This  is  readily  accomplished  in  the  algorithm  by  special  initialization  of  the  wave  function. 
In  particular,  given  the  complex  BEC  ground  state  wave  function^ (x,r  =  0)  =  a(x,0)  +  j3(x,0) , 
we  simply  choose 

a(x,0)=  ilmij/(x,0)  ,  /?(x,0)=  Rey/(x,0)  (25) 

We  also  need  to  generalize  the  phase  rotation  operator  Q.  to  the  non-diagonal  potential  matrix 


a  = 


cos(|i/a|2) 

-/sin(yx2) 

-/sin^l2) 

cos(|v^|2) 

so  that 


£(6  •  W))y  =exp[-i|v/" 


¥ 


(26) 


(27) 


on  setting  e  1 . 


5.  Scaling  of  the  Qubit  Lattice  Algorithm  on  Various  Architectures 

While  all  our  production  runs  have  been  performed  at  the  DoD  supercomputer  sites,  the 
scaling  experiments  where  performed  on  the  larger  DoE  supercomputers  at  Argonne  (IBM/Blue 
Gene),  Oak  Ridge  (the  CRAY  Jaguarpf  Titan)  and  Oakland  (NERSC  CRAY  Hopper).  In 
particular,  using  a  coupled  BEC  qubit  unitary  algorithm  we  have  achieved  excellent  strong  and 
weak  scaling. 

In  strong  scaling,  the  grid  is  held  fixed  while  one  increasing  the  number  of  processors. 
Ideal  scaling  would  result  in  a  proportional  decreasing  in  the  wall  clock  time.  On  the  Oak  Ridge 
CRAY,  we  have  close  to  perfect  scaling  to  276  480  cores  tested. 


Strong  Scaling:  Grid  84003,CRAYX 

K-6  (Titan/Jaguarpf) 

#processors 

Wall-clock  (s) 

Speed-up 

Ideal  Speed-up 

96  000 

1422.2 

1 

1 

120  000 

1087.9 

1.25 

1.25 

150  000 

872.9 

1.63 

1.56 

180  000 

729.3 

1.95 

1.88 

210  000 

643.5 

2.21 

2.19 

216  000 

628.6 

2.26 

2.25 
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Strong  Scaling:  Grid  96003,CRAYX 

K-6  (Titan/Jaguarpf) 

#processors 

Wall-clock  (s) 

Speed-up 

Ideal  Speed-up 

120  000 

1674.2 

1 

1 

150  000 

1337.4 

1.25 

1.25 

240  000 

834.4 

2.01 

2.00 

262  144 

769.2 

2.18 

2.18 

276  480 

739.4 

2.26 

2.30 

Fig.  1  Strong  scaling  of  the  coupled-BEC  code  on  Jaguarpf/Titan  CRAY  XK-6.  The  red  dots  are 
ideal  scaling  while  the  blue  dots  are  our  qubit  unitary  lattice  code. 
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On  the  IBM/Blue  Gene  Intrepid,  we  find  similar  excellent  strong  scaling  -  tested  up  to  163  840 
cores: 


Strong  Scaling:  Grid  48003 ,  IBM/Blue  Gene  (Intrepid) 


#processors 

Wall-clock  (s) 

Speed-up 

Ideal  Speed-up 

65  536 

1134.5 

1 

1 

98  304 

765.0 

1.48 

1.5 

131  072 

561.3 

2.02 

2.0 

163  840 

455.0 

2.49 

2.5 

Weak  scaling  is  when  one  changes  the  grid  proportionally  to  the  number  of  cores  so  that 
each  processors  does  the  same  number  of  computations.  Ideal  scaling  would  result  in  the  wall 
clock  being  constant.  Again  our  codes  scale  excellent  under  weak  scaling  on  both  the  CRAY  XK-6 
and  the  IBM/Blue  Gene 


Weak  Scaling:  CRAY  XK-6  (Titan/Jaguarpf) 


#processors 

Grid 

Wall  Clock  (s) 

Ideal  (s) 

1  728 

1200 3 

223.5 

223.5 

13  824 

2400 3 

224.8 

223.5 

46  656 

3 600 3 

234.0 

223.5 

110  592 

4800 3 

225.7 

223.5 

216  000 

6000 3 

231.6 

223.5 

Weak  Scaling  :  IBM/Blue  Gene  (Intrepid) 


#processors 

Grid 

Wall  Clock  (s) 

Ideal  (s) 

64 

400 3 

425.9 

425.9 

512 

800 3 

426.9 

425.9 

4  096 

1600 3 

425.5 

425.9 

32  768 

3200 3 

421.3 

425.9 

Work  has  also  been  performed  on  restart  i/o  using  Adios  on  the  CRAY  XK-6,  and  again 
we  are  seeing  extremely  high  GB/s  being  achieved: 


Restart  i/o  (Adios)  on  the  CRAY  XK-6  (Titan/Jaguarpf) 


Grid 

#processors 

Write  (s) 

GB/s  (write) 

Read  (s) 

GB/s  (read) 

1200 3 

1  728 

3.9 

13.2 

2.5 

20.6 

2400 3 

13  824 

18.0 

22.9 

22.3 

18.5 

3 600 3 

46  656 

35.0 

39.7 

36.0 

38.6 

4800 3 

110  592 

75.8 

43.5 

74.9 

44.0 
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6.  3D  Quantum  Turbulence 

Following  Berloff,  we  can  determine  the  Pade  approximant  steady  state  quantum  line 
vortex  solution  to  the  GP  Eq.  (24).  In  particular  for  winding  number  n  =  1 


1 


1 1  ar‘ 


(l2  +  ar2 ) 


g[384  +  ar2(l28  +  llar2) 
— »  g_1/2  as  r  — >  oo  ,  and  k 


(28) 


— >  rall2g  1/2  as  r  — >  0  . 


where  asymptotically: 

Since  the  GP  Eq.  (1 1)  is  a  Hamiltonian  system,  the  total  energy  of  the  BEC  is  a  constant: 

E  =  Ekin  +  Equ  +  Emt  =  COnst  ■  (29) 


(30) 


(31) 


where  the  total  kinetic  energy 

E kin  —  2  J"  ^  ^  I  ^  I  ’ 

the  quantum  energy  (related  to  the  unusual  quantum  stress  tensor) 

Equ=2jd3x\V^f  , 
and  the  interaction  energy 

Eint  =  as\ d3x  P 2  ■  (32) 

The  kinetic  energy  can  be  further  subdivided  into  the  compressible  and  incompressible 
components  using  the  Helmholtz  decomposition: 

rp  jpCOMP  .  rp  INCOMP  /qq\ 

Ekin  ~  Ekin  +  Ekin  \66) 

However,  since  the  quantum  vortex  is  a  topological  singularity  we  need  to  regularize  the  standard 
definition  of  velocity  and  vorticity  fields.  In  particular  we  work  with  the  density  weighted  velocity 

field  s[pns  to  regularize  the  superfluid  velocity  at  the  quantum  vortex  core. 

The  development  of  quantum  turbulence  can  be  readily  seen  in  the  time  evolution  of  the 
quantum  vortex  cores  in  a  run  on  a  40323  grid,  Fig.  8 


(a)  t=  0  (b)  t  =  400  (c)  t  =  57  500 

Fig.  8  A  zoomed-in  section  of  the  GP  time  evolution  of  the  isosurfaces  of  2  quantum  line  vortex 

cores  (from  a  total  of  12  such  vortices)  on  a  spatial  grid  of  40323 .  (a)  Att  =  0 ,  the  vortex  core 
isosurfaces  for  initial  winding  number  6  -  which  is  very  unstable.  By  t  =  400,  (b),  the  cores  have 
expanded  and  fragmented  into  n  =  1  vortices,  (c)  Vortex  tangle  is  evident  at  time  t  =  57  500. 
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Of  much  interest  is  the  energy  spectra  for  3D  quantum  turbulence.  Because  of  the 
excellent  parallelization  of  our  qubit  algorithm,  we  have  received  over  6  CAP  special  user 
allocations  over  these  past  7  years.  This  has  permitted  us  to  perform  production  runs  on  grids  not 
achievable  by  standard  direct  CFD  solvers.  In  particular,  we  have  examined  the  energy  spectrum 
on  a  57603  grid.  Because  ofcomputational  time  constraints,  the  spectral  analysis  was  performed 
on  MRAP  over  a  2  week  period  just  after  it  was  acquired  and  just  before  it  became  a  classified 
machine.  Hence  we  could  only  analyze  the  total  kinetic  energy  spectrum,  Fig.  9. 


Total  Kinetic  Energy  EK1N  grid  57603 


Fig.  9  The  total  kinetic  energy  spectrum  Ekin(k )  from  a  57603  grid  simulation.  There  are  3  clear 
power  law  spectral  exponents:  For  large  k  (800  <k<  2500)  Ekin  (k)  ~  k~3 :  for  intermediate  k  (300 
<k<  700)  a  steep  Ekin  (k)  »  r7'8 ;  while  for  small  k  (100  <k<  200)  Ekin  (k)  »  k~5'3 . 


While  the  spectral  exponents  for  the  medium  and  large-k  regimes  are  very  robust,  there  are 
considerable  fluctuations  in  the  small-k  spectral  exponent  (an  exponent  of  -1.30  if  one  considers 
the  region  30  <  k  <  175). 

A  detailed  look  at  the  individual  energy  spectra  was  performed  on  grids  of  3072  i  and 
the  effect  of  initial  winding  number  on  the  line  vortices,  Figs.  10-11.  At  t  =  0,  over  99%  of  the 
kinetic  energy  is  in  the  incompressible  mode.  However,  this  rapidly  changes,  and  even  by  t  =  10  K, 

the  triple  cascade  spectral  ranges  are  readily  seen  in  both  Ek°MP  (&)  and  Equ  (&) .  However,  the 

EIkfnCOMP(k)  spectrum  shows  only  a  k~ 3  spectrum  over  the  whole  k-range:  30  <  k  <  1500.  As 

expected,  the  n  =  2  spectra  are  less  noisy  for  the  initial  vortices  with  winding  number  n  =  2  (Fig. 
11b  compared  to  Fig.  11a). 
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Fig.  10  initial  energy  spectra  Fig.  11  energy  spectra  at  t  =  10  000 

(a)  non-degenerate  vortices  with  winding  number  n=1,  (b)  degenerate  vortices  with  n  =  2. 

Color  scheme:  blue  -  eZcomp  (k) ,  red  -  e£°mp  (k)  ,  gold-  Equ(k). 

Using  linear  regression  in  fixed  wave  number  k  windows,  we  have  computed  the  time  averaged 
energy  spectral  exponents  (a) ,  and  their  average  deviation  ±(<j)  .  The  time  average  is 

performed  on  60  data  points  from  t  =  24  000  to  t  =  48  000  in  steps  of  At  =  400 

Incompressible  Kinetic  Energy  Spectrum 


k-band 

(a) :  E{k)~k~^ 

±((7> 

15  <  k  <  90 

2.72 

0.21 

180  <  k  <  280 

3.377 

0.021 

350  <  k  <  900 

3.006 

0.011 
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Quantum  Energy  Spectrum 


k-band 

(a) :  E{k)~k~^ 

±((7> 

15  <  k  <  90 

1.70 

0.17 

180  <  k  <  280 

7.802 

0.058 

350  <  k  <  900 

3.033 

0.018 

Total  Kinetic  Energy  S 

pectrum 

k-band 

(a) :  E[k)~k~^ 

±((7> 

15  <  k  <  90 

1.66 

0.17 

180  <  k  <  280 

8.527 

0.064 

350  <  k  <  900 

3.042 

0.020 

We  conjecture  that  for  large  k  >  350  one  is  basically  looking  at  the  intravortex  spectrum  of  a 

“single”  vortex  as  the  incompressible  kinetic  energy  spectrum  is  a  very  robust  k~ 3 .  However,  what 
makes  it  more  complicated  is  that  both  the  compressible  and  quantum  energy  spectra  themselves 

also  exhibit  a  strong  k~ 3  -  a  spectral  property  clearly  not  exhibited  by  a  single  linear  quantum 
vortex  (c.f.,  Fig.  10).  Around  £  =  350,  both  the  compressible  and  quantum  energy  spectra  lift 
sharply  away  from  the  k~ 3  spectrum  as  we  move  to  smaller  k.  Around  this  sharp  rise,  the 
incompressible  kinetic  energy  spectrum  moves  significantly  away  from  k~ 3  to  k~ 3377  for  k  <  280. 
Whether  this  higher  exponent  is  due  to  Kelvin  waves  on  the  vortices  needs  further  investigation. 
The  small  k  spectral  region  is  also  of  great  interest,  since  it  is  tempting  to  identify  this  region  with 

the  classical  Kolmogorov  A:_5/3  kinetic  energy  spectrum  of  fluid  turbulence.  In  our  quantum 
turbulence  simulations,  the  incompressible  kinetic  energy  spectrum  does  not  exhibit  a  k~5'3  for 
small  k  -  but  the  total  kinetic  energy  does.  Interestingly,  in  subgrid  closure  models  of  compressible 
classical  fluid  turbulence,  it  is  the  total  subgrid  kinetic  energy  spectrum  that  is  assumed  to  satisfy 
the  Kolmogorov  k~5'3  spectrum. 


7.  2D  Quantum  Turbulence 

There  has  been  some  interest  to  see  whether  2D  quantum  turbulence  would  show  energy 
cascade  properties  like  in  3D  quantum  turbulence  or  whether  it  would  follow  classical  physics:  in 
incompressible  classical  fluid  turbulence,  there  is  direct  energy  transfer  from  low-k  to  high-k  modes 
for  3D  turbulence  while  there  is  an  inverse  cascade  of  energy  from  high-k  to  low-k  modes  as  in  2D 
turbulence.  Because  the  quantum  vortex  is  singular  at  the  core,  we  renormalize  the  vorticity  by 
weighting  it  with  the  density 


co, 


=  [Vx(V7“ 


(34) 


Thus  coq  exhibits  the  effects  of  density  variations  ( V^/p  xu,  •  z )  as  well  as  the  2D  location  of  the 
vorticies  ^/pV  xusx-z  =  ^/p  r<5(x-x0)  .  We  also  consider  the  density-weighted  enstrophy 


z  = 


(35) 
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It  is  interesting  how  vortices  are  induced  in  BECs.  Experimentally,  with  sufficient  rotation 
of  the  BEC,  a  lattice  of  vortices  is  induced  within  the  BEC.  Vortices  can  also  be  induced  by  strong 
Kelvin-Helmholtz  instabilities,  as  we  have  seen  from  our  simulations,  Fig.  12. 


# 


\ 


(t>)  0 .  t=l(M>00 


(n)  0 ,  t as  10500 


Fig.  12  The  generation  of  vortex-antivortex  pairs  in  the  quantum  Helmholtz  instability  in  2D.  The 
upper  plots  are  for  the  vorticity  coq  (x,t)  while  the  lower  plots  are  the  corresponding  phases  of  the 

wave  function,  (a)  at  t  =  10500,  there  are  4  vortices  present  in  this  2D  spatial  region:  in  the  upper 
plot,  the  red  localized  points  are  for  vortices  rotating  anticlockwise,  while  the  blue  localized  points 
are  for  vortices  rotating  clockwise.  In  the  lower  plots  of  the  wave  function  phases,  blue 
corresponds  to  phase  0  while  red  corresponds  to  phase  In .  The  location  of  the  branch  point 
singularities  and  the  branch  cuts  are  very  evident  (winding  number  1),  corresponding  to  the  strong 
dots  seen  in  the  top  frame.  There  is  strong  shear  between  the  two  like-wise  rotating  vortices, 
leading  to  quantum  Kelvin-Helmholtz  instability  and  the  generation  of  a  pair  of  clockwise-and- 
anticlockwise  point  vortices,  as  pointed  out  by  the  arrows  in  (b).  While  the  circulation  of  each 
vortex  is  the  same,  the  density  weighting  around  these  vortices  leads  to  different  strengths  of  the 
density-weighted  vorticity  coq . 


Another  interesting  simulation  considers  the  evolution  of  2D  quantum  turbulence  from  the 
initial  condition  of  constant  density  p- const.,  but  random  phase.  Thus  initially  there  are  no 
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quantum  vortices  (since  at  a  quantum  vortex  core  p  =  0).  Since  the  superfluid  velocity  is  the 
gradient  of  this  random  phase,  the  initial  state  is  very  unstable,  and  so  quantum  vortices  will  rapidly 
develop. 


(d)  1  =  20800 


(e)  t  =  21000=  Tp/2 


(f)  1  =  21100 
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Fig.  13  The  time  evolution  of  the  density-weighted  vorticity  |oj(/(x,/)|  from  a  constant  density, 

random  phase  initial  condition.  At  a  quantum  vortex  core  a>q~8(x-xi).  (a)  t=0,  no  quantum 

vorticies  since  p(x,t  =  0)^0.  (b)-(d)  quantum  vortices  (  8  -  function  spikes  )  and 

compressional  waves,  (e)  the  onset  of  a  semi-Poincare  recurrence  in  which  the  vorticity  tends  to 
its  initial  state.  There  are  no  quantum  vortices,  but  there  are  compressional  waves,  (f)  The 
reappearance  of  quantum  vortices. 


Because  the  GP  system,  Eq.  (11),  is  a  Hamiltonian  system  there  must  be  a  Poincare 
recurrence  of  the  initial  conditions.  However,  for  almost  all  continuous  systems,  this  Poincare 
recurrence  time  effectively  has  Tp  — » °° .  For  BECS  under  certain  conditions  -  particularly  when 

the  internal  energy  «  total  kinetic  energy  -  we  find  remarkably  short  Poincare  recurrence  times. 
Moreover,  we  find  point  inversion  symmetry  of  the  initial  conditions  at  Tp/2.  This  is  reminiscent 
of  the  2D  Arnold  Cat  Map  on  the  [0,1]  square: 

r2x  +  y^\ 

— >  ,mod  1  (36) 


yyj  yy  j 

The  Arnold  Cat  Map  is  invertible,  area  preserving,  mixing  and  ergodic  with  relative  short  Poincare 
recurrence  time  Tp .  Indeed  for  pixel  resolution  ( 300  x  300 ) 


(a)  t=0  (b)  t=150  (c)  t=300 


Fig.  14  Poincare  recurrence  of  initial  conditions  for  a  particular  pixel  resolution,  300  x  300 ,  for  the 
Arnold  Cat  Map. 


Yet  for  other  pixel  resolutions,  one  sees  a  point  inversion  symmetry  at  TP  /2 : 


(a)  t=0  (b)  t=57  (c)  t=l  14 

Fig.  15  Poincare  recurrence  of  initial  conditions  for  a  pixel  resolution,  74x74,  for  the  Arnold  Cat 
Map.  Point  symmetry  inversion  of  the  initial  conditions  are  now  found  after  Tp  /  2  iterations. 
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The  precursor  to  the  occurrence  of  a  short  Poincare  recurrence  time  for  the  GP  system,  Eq.  (11), 
can  be  seen  in  the  time  evolution  of  the  kinetic,  quantum  and  total  energies  as  well  as  the 
enstrophy: 


(a)  kinetic,  quantum  and  total  energy  (b)  enstrophy  ,  Eq.  (35). 

Fig.  16  The  spikes  in  the  time  evolution  of  the  (a)  kinetic  and  quantum  energies,  and  (b)  enstrophy 
indicate  the  times  of  semi-Poincare  recurrence  for  random  phase,  constant  density  initial 
conditions.  The  strong  dips  in  the  enstrophy  indicate  the  depletion  of  quantum  vortices. 


However,  for  stronger  internal  energies,  where 


one  finds  the  destruction  of  the  short  Poincare  recurrence  time,  Fig.  17.  The  average  (  )  is  a  time 
average. 


(c)  Energy  plot  for  7  =  1.33  x  10_1  (d)  Enstrophy  plot  for  7  =  1.33  x  10-1 
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Fig.  17  The  loss  of  the  short-Poincare  recurrence  time  as  the  initial  internal  energy  is  increased. 

The  spectral  properties  of  2D  quantum  turbulence  is  complex.  However,  all  indications 
are  that  we  do  not  see  inverse  cascade  of  energy  to  small  scales.  If  we  concentrate  on  the 
incompressible  kinetic  energy  spectrum,  then  we  find  the  following  spectral  exponents  for  2  runs: 
one  with  initial  vortices  with  winding  number  1,  Fig.  18a,  and  the  other  with  winding  number  2,  Fig. 
18b. 


Plot  of  sIC  (t)  for  (a)  winding  number  1  (b)  winding  number  2 

Fig.  18  Incompressible  kinetic  energy  spectrum  exponent  E^OMP  ( k,t)~  ks,c(t) .  The  red  line  is 

for  the  incompressible  spectrum  with  sIC  =  -3 :  E^OMP  (k)  ~  k~3 .  The  steep  increase  in  the 

spectral  exponent  can  be  traced  to  the  strong  depletion  of  quantum  vortices  in  the  system  during 
that  period  of  evolution  of  the  quantum  turbulence.  For  the  more  turbulent  system  with  winding 
number  2,  one  does  not  see  such  a  systematic  destruction  of  vortices  and  hence  the 
incompressible  kinetic  energy  spectrum  has  only  small  fluctuations  about  sIC  =  -3 ,  (b). 


phases  (p(x,y-,t):  (al)  t=  24400 


(bl)  t  =  24500 


(cl)  t=  24600 
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Fig.  19  The  evolution  of  the  incompressible  kinetic  energy  spectrum  for  two  neighboring  times 
about  the  peaking  in  the  exponent  E^OMP  (k)  ~  k~5882  at  t  =  24500  ,  with  the  strong  kink  (b) 
emphasized  by  the  circle  in  the  large-k  region.  The  corresponding  blow-up  of  the  phases 
(p(x,y;t)in  a  very  limited  spatial  subdomain  where  there  are  6  quantum  vortices  at  t  =  24400, 
(al),  NO  vortices  at  t  =  24500,  (bl),  and  the  re-emergence  of  4  quantum  vortices  at  t  =  24600 . 

A  closer  look  at  what  is  happening  about  these  sharp  spectral  exponent  increases  in  the 
incompressible  kinetic  energy  is  given  in  Fig.  19.  At  t  =  24400,  (a)  and  (al),  the  incompressible 

kinetic  energy  spectrum  shows  a  long  range  of  E^OMP  (k)  ~  k~3 ,  and  we  note  in  the  spatial 
subdomain  that  there  are  6  vortices/branch  points  with  their  corresponding  branch  cuts.  By  t  = 

24500,  (b)  and  (bl),  the  incompressible  spectrum  has  steeped  considerably  to  E^OMP(k)  ~  k~6 
with  a  kink  for  large-k  just  before  the  spectrum  moves  away  from  a  power  law.  An  examination  of 
the  same  spatial  subdomain  as  in  (al)  indicates  that  the  quantum  vortices  have  been  destroyed  as 
they  are  no  branch  points  in  (bl).  These  quantum  vortices  re-emerge  at  t  -  24600  (actually  now  4 

quantum  vortices)  and  the  incompressible  spectrum  relaxes  again  to  E^OMP  (k)  ~  k~3A . 

The  2D  quantum  spectrum  has  also  been  found,  for  strong  winding  number 
degeneracies,  to  resemble  the  3D  quantum  spectrum  -  with  the  compressible  kinetic  energy 
spectrum  exhibiting  a  triple  cascade  region,  Fig.  20 
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(a)  compressible  kinetic  energy  spectrum  (b)  incompressible  kinetic  energy  spectrum. 
Fig.  20  The  kinetic  energy  spectra  for  an  initial  12  quantum  vortex  state  of  degeneracy  6  on  a  grid 

of  81922.  (a)  The  compressible  spectrum  exhibits  a  triple  cascade  behavior  with  exponents 
sc  =  -1 .83  for  small  k,  sc  =  -8 . 1  for  intermediate  k  and  sc  =  - 2.85  for  large  k.  (b)  The 
incompressible  spectrum  has  a  sIC  =-2.93  exponent  (red),  while  for  an  intermediate  range  we 
see  sIC  =  -4.00  (green). 


8.  Skyrmions 

As  one  moves  away  from  scalar  topological  defects  to  higher  dimensional  ones,  we  have 
started  to  focus  on  skyrmions.  These  are  localized  in  space  with  quantized  topological  charge  - 
not  unlike  solitons  -  and  which  can  undergo  interactions  and  form  ordered  phases  and  undergo 
phase  transitions.  Such  vortices  have  been  produced  experimentally  using  the  F  =  1  Na 
condensate.  Skyrme  (1950’s)  initially  introduced  his  ideas  to  explain  neutrons  and  protons  as 
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discrete  entities  emerging  from  a  continuous  field.  When  quarks  entered  into  the  vocabulary  in  the 
1960’s,  skrymions  were  neglected  and  then  were  again  overlooked  when  string  theory  became  the 
vogue.  However,  with  the  discovery  of  the  quantum  Hall  effect  in  the  1990’s  interest  has  been 
rekindled  in  skyrmions  because  of  the  topological  properties  in  the  Hall  effect  and  in  skyrmions. 
We  consider  the  skyrmion  field  as  a  spinor  (in  cylindrical  coordinates) 
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The  vortex  cores  as  well  as  the  asymptotic  wave  function  isosurfaces  are  shown  in  Fig.  21 


(a)  vortex  core  isosurfaces  \\j/+\2  =  0.04  =  \x//_\2 ,  (b)  asymptotic  isosurfaces  |y/+|2  =  0.9  =  \i//_\2 
Fig.  21  Isosurfaces  of  the  2-component  spinor  skyrmion  field,  (a)  The  vortex  cores  are  a  ring 
vortex  for  y/+  (green) ,  and  for  y/_  (red)  a  linear  vortex  along  the  vertical  axis  which  then  closes 
asymptotically  over  this  apple  core  at  infinity,  (b)  The  asymptotic  surfaces  interchange  roles. 


Note  that 
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,  as  x  — >  °° 


The  skyrmion  evolution  is  given  by  the  solution  of  the  spinor  GP  equation 

i^=  -  V!V  +a(S{|v/+|2-|^_|2}-l)T  (38) 

Since  we  are  interested  in  quantum  turbulence,  we  rescale  the  skyrmion  initial  condition  ,  Eq.  (37), 
so  that  it  is  not  a  quasi-eigenfunction  for  the  spinor  GP  Eq.  (38).  Hence  we  see  the  quantum 
vortex  core  of  the  spinor  expand,  then  interact  with  the  periodic  walls  to  form  a  percolating  lattice  of 
many  skyrmions,  Fig.  22. 
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(c)  t  =  1000  (d)  t  =  2500 


(e)  1=  7000  (f)  1=  10  500 
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(k)  t  =  28  000  (l)  t  =  29  000 
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Fig.  22  The  evolution  of  the  spinor  quantum  core  of  a  skyrmion  under  periodic  boundary 
conditions.  The  single  skyrmion  rapidly  breaks  into  a  lattice  of  percolating  skyrmions.  These 
isosurfaces  are  formed  by  culling  the  front  interface  to  expose  the  innards.  The  initial  vortex  core  is 
a  ring  vortex  for  y/+  (green) ,  and  for  yr  (red)  a  linear  vortex  along  the  vertical  axis  which  then 
closes  asymptotically  over  this  apple  core. 


A  snapshot  of  the  spectra  for  the  linear-toroidal  vortex  core  y/_  (red  in  Fig.  22)  are  shown  in  Fig.  23 

10*°r - - - - - - - - - j  10Wr - - - - - - , 


(a)  t=  0 


(b)  t=  IK 


(c)  t=  24000  (d)  t=  27000 

Fig.  23  Kinetic  energy  spectrum  (blue),  quantum  energy  spectrum  (red)  and  the  internal  energy 
spectrum  (green)  for  the  isosurfaces  of  the  spinor  component  \y/_\ .  The  dashed  line  in  (b)  is  the 

k~ 3  spectrum. 


The  spectrum  for  the  vortex  ring  core  y/  (green  in  Fig.  22)  is  different,  and  shown  in  Fig.  24 


(a)  t=  0 


(b)  kinetic  energy  spectrum  at  3  times 
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k  k 


(c)  t=  14000  (d)  t=  27000 

Fig.  24  Kinetic  energy  spectrum  (blue),  quantum  energy  spectrum  (red)  and  the  internal  energy 
spectrum  (green)  for  the  isosurfaces  of  the  spinor  component  \\j/+\ .  While  the  initial  spectra  for 

kinetic  and  quantum  energies  exhibit  the  ubiquitous  k~ 3  spectrum,  this  is  rapidly  destroyed  and 
followed  by  steep  non-power  law  dacays.  (b)  exhibits  the  rapid  destruction  of  this  k~ 3  spectrum 
for  the  kinetic  energy:  t  =  0  (blue),  t=  1000  (red) ,  t  =  2000  (green). 
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